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Abstract 

We consider pricing weather derivatives for use as protection against weather extremes. The 
method described utihzes results from spatial statistics and extreme value theory to first 
model extremes in the weather as a max-stable process, and then use these models to sim- 
ulate payments for a general collection of weather derivatives. These simulations capture 
the spatial dependence of payments. Incorporating results from catastrophe ratemaking, we 
show how this method can be used to compute risk loads and premiums for weather deriva- 
tives which are renewal- additive. 

Keywords: extreme value, generalized extreme value distribution, max-stable process, 
renewal-additive, weather derivative 

1. Introduction 

Weather derivatives are contingent contracts whose payments are determined by the 
difference between some underlying weather measurement and a pre-specified strike value. 
They provide a useful risk management tool for any party facing weather risk. They also 
provide investments which are often uncorrelated with more traditional financial instruments, 
allowing investors to diversify. The first weather derivative was developed in 1996, and by 
1999 derivatives and their options were being traded on the Chicago Mercantile Exchange 
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(Kunreauther and Michel- Kerj an, 2009). 

Richards et. al. (2004) give a hst of 5 elements common to all weather derivatives. These 
include (a) an underlying weather index, (b) a well-defined time period, (c) the weather 
station used for reporting, (d) the payment attached to the index vahic, and (e) the strike 
value which first triggers payment. The intention is for the buyer of the derivative to be 
compensated by the seller for amounts which roughly correspond to actual business losses. 
Ideally, these losses are perfectly correlated with the payments of the weather derivative, 
though in practice this is rarely achieved. Tailoring the contract to the specific needs of one 
buyer reduces its general appeal in a secondary market, and thus lowers the value of the 
contract. 

Weather derivatives offer benefits to the buyer and seller not found in traditional insur- 
ance. The buyer docs not need to have an insurable interest, nor do they need to demonstrate 
an actual loss to receive payment. The loss payment itself is generally proportional to the 
difference between the weather index and strike value. Furthermore, weather derivatives of- 
fer the buyer the opportunity to sell (or even buy back) the contract in a secondary market 
such as the Chicago Mercantile Exchange, or over the counter. There are no such markets 
for traditional insurance products. Derivatives are also in general more lightly regulated, 
and payments are often considered taxable. 

Prom the seller's point of view, weather derivatives offer several advantages over tradi- 
tional insurance. Weather derivatives avoid the higher administrative and loss adjustment 
expenses of insurance contracts. They also eliminate concern for moral hazard, morale haz- 
ard, and fraud, as the event triggering the payment is easily verified and completely beyond 
the buyer's control. When used to insure crops, weather derivatives help reduce the per- 
ceived information asymmetry associated with crop insurance, wherein farmers often have 
more information about their individual risk than insurers (Goodwin and Smith, 1995). 

In this paper we consider the use of weather derivatives to provide protection against 
high impact, low probability business losses caused by extremes in the weather. Though 
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not technically insurance, we model losses and price derivatives using actuarial techniques 
originally developed to price insurance. The weather derivatives we consider define some 
payment L = L{M] s, t), where M is the unknown weather random variable, and s and t are 

the pre-specified strike and limit values (occasionally we only write L or L{m) to refer to 
the loss to simplify notation). Examples of three types of derivatives with payments based 
on high exceedances include 

1. L = q; if {M > s} and otherwise 

2. L — /3 ■ {m — s) a {M > s} and otherwise 

3. L = (5 ■ {m — s) when {s < M <t} and L = {3 -{t — s) when {M > t}, and otherwise, 

where a and j3 are dollar values, and m is the realization of random variable M. The first 
provides a fiat payment whenever the event {M > s} occurs, the second provides a propor- 
tional payment based on the difference (m — s), while the third limits the total payment. 
Unlike most derivatives, weather derivatives do not have an underlying tradable asset, and 
thus many pricing approaches based on financial theory are inappropriate. Jewson and Brix 
(2005) provide an excellent reference and discussion of pricing techniques for weather deriva- 
tives. The pricing approach we take is based on computing expected losses and expected 
loss variability. The first step is to compute the expected payout E{L) — J L{m)g{m) dm, 
where g{m) is the density function of weather variable M. For the three derivatives, shown, 
expected payments are 

1. E{L) ^ J^a- g{m) dm ^ a ■ P{M > s) 

2. E{L) = /3 ■ (m — s)g{m) dm 

3. E{L) = ^ • (m - s)g{m) dm + p-(t-s)- P{M > t). 

When viewed from the point of view of insurance, the quantities E{L) arc the pure 
premiums of the contracts. In the absence of any expenses, profit, risk loadings and time- 
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value financial considerations, this is the price of the contract for the buyer. Such contracts 
are fairly straightforward to price once one has an accurate estimate of the density function 
g{m) in the region where m> s. Further, one can compute all second moments as 

1. E{L'') = /~ . g{m) dm^a^- P{M > s) 

2. E{L^) = J^{/3 ■ (m - s)yg{m) dm 

3. E{L^) = f^{/3 ■ (m - s)}^g{m) dm + {/3 ■ {t - s)}^ • P(M > t), 

and from these compute the variance var(L) = E{LF') — {E{L)y. This information is often 
incorporated into the premium by adding a risk load R{L), which is added to the pure 
premium as P = E{L) + R{L). Risk loads are meant to account for the additional risk 
taken on by writing derivatives with larger variability of losses. Commonly, risk loads are 
a function of the variance (or standard deviation) of loss (Feldblum, 1990). Mango (1998) 
describes several common risk loadings such as R{L) = X ■ ^var(L), or R{L) = A • var(L), 
where A is a dollar amount chosen to satisfy some risk tolerance criteria. 

Next, consider a portfolio of K weather derivatives, with aggregate payment L — Li + 
... + Lk. The expected aggregate payment is simply the sum of the individual expected 
payments, 

K 

E{L) = Y,E{Lk). 

k=l 

However, when we allow for possible dependence among contracts, the variance of the ag- 
gregate payment is 

K K-l K 

var(L) = ^ var(Lfe) + XI ^ ' ^^^i^k, Lk') (1) 

fc=l k=l k'=k+l 

The first issue when pricing weather derivatives (extreme or otherwise) is to properly ad- 
dress the positive correlation among contracts, with its resulting impact on aggregate loss 
variability as shown in equation 1. It is easy to envision positively correlated payments in 
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a portfolio of weather derivatives, since weather variables are often positively correlated in 
space. This correlation implies the variance of the aggregate payment exceeds the sum of 
individual payment variances, and thus risk loadings priced individually would be insufficient 
for the portfolio whole. 

Next, consider the challenges when focusing on extreme weather events. Examples of 
such events may include the maximum daily temperature exceeding some high threshold, 
the minimum daily temperature falling below some low threshold, or the minimum monthly 
rainfall falling below some low threshold. These events can often be written as {maxF > s} 
or {miny < s} where Y is some weather- related random variable and s is a pre-specified 
strike value. In all cases we are defining a contract not based on "typical" weather patterns 
of temperature or precipitation, but on extremes. What is needed to accurately price these 
contracts is the distribution function of extreme events. Furthermore, when one consid- 
ers a collection of K derivatives defined at different locations, one must carefully consider 
if the dependence of extremes is different from the dependence exhibited by non-extreme 
events. Putting the two issues together, what is needed to price weather derivatives for 
extreme events is a model that (1) directly targets extremes, and (2) properly incorporates 
the spatial correlation of weather extremes. One could further extend this to models which 
incorporate time dependence of extremes as well. However, in this paper we focus on the 
spatial dependence of weather derivatives for extremes, as that is the largest omission of 
current methodology. 

We begin with the problem of pricing a single weather derivative in the first section. This 
is handled through the Generalized Extreme Value distribution (Embrechts, et. al., 1999), 
which is the only permissible limit of the maxima of independent, identically distributed 
univariate random variables. We demonstrate the approach by pricing a weather derivative 
based on extreme summer temperatures in Phoenix, Arizona. In the third section, we intro- 
duce models used in spatial statistics and spatial extremes, which will serve as a foundation 
for the fourth and fifth sections. There, we extend our weather derivative pricing model 
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to multiple locations through the use of max-stable processes. These processes capture de- 
pendence in spatial extremes. Through large numbers of simulations, we can estimate all 
marginal variances, covariances, and other quantities of interest when pricing a portfolio of 
weather derivatives. This information is ultimately incorporated into risk loads added to the 
pure premiums. From the method presented, pure premiums and risk loads for a collection 
of K spatially dependent weather derivatives can be obtained. An application is shown in 
section 6. 



2. Pricing a Weather Derivative for an Extreme Event at a Single Location 

2.1. The Generalized Extreme Value Distribution 

We begin by introducing the model used for modeling maxima (minima can always be 
rewritten as min(Yi, y„) = — max(— Yi, —1^), so there is no loss in generality when one 
considers only maxima). Let l^i, be independent and identically distributed univariate 
random variables with some distribution function F, and let M„ = max(Yi, be the 

maximum. If M„ converges to a non-degenerate distribution under re-normalization as 

^" < m ) = F'*(a„m -|- — >■ G{m) as n — >■ oo 



On / 

for some sequences a„ and 6„, then G must be a member of the Generalized Extreme Value 
(GEV) family, with distribution function 

G(m)=exp|-(^l + e^) '^'|. (2) 

Here a+ = max(a, 0), and cr, and ^ are the location, scale, and shape parameters, respec- 
tively (Coles, 2001). The sign of the shape parameter ^ corresponds to the three classical 
extreme value distributions: ^ > is Prechet, { < is WeibuU, and ^ ^ (in the hmit) is 
Gumbel. The Frechet case corresponds to a heavy tailed distribution, Gumbel is intermedi- 
ate, and WeibuU has a bounded upper limit. 
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In practice, it isn't necessary to worry about specifying sequences 6„ and a„ due to the 
property of max-stability: li Yi, ...,Yn are independent and identically distributed from G, 
then max(Y'i, ¥„) also has the same distribution with only a change in location and scale, 
as C^ly) = G{Any + Bn) for constants An and -B„. A distribution is a member of the 
Generalized Extreme Value family if and only if it is max-stable (Lcadbcttcr ct. al., 1983). 
A practical consequence of this property is that if one changes the block size (from monthly 
maxima to annual maxima, for instance) and fits a new model, new estimates for the three 
GEV parameters (//, a, ^) are obtained, but the model is still GEV. 

A special case of the Generahzed Extreme Value family is the unit-Prechet, with distri- 
bution function G{m) — exp(— 1/m). Any member of the Generalized Extreme Value family 
may be transformed to have unit-Prechet margins as follows: if Y has a Generalized Extreme 
Value distribution with range < F < oo, then a new variable U may be defined as 



and U has unit-Prechet margins. If the parameters are unknown, they may first be estimated 
and then the transformation to U is taken. When we model multivariate or spatial extremes, 
there is no loss in generality when one assumes the margins are all unit-Frechct. In practice, 
one would first estimate all marginal distributions and transform to unit-Prechet, then in a 
second step analyze the spatial dependence. 

The GEV model can be fit to observed data using maximum hkelihood estimation. Call 
the parameter vector 0. This parameter can be as simple as three fixed parameters, as 
(f) — {ii,a,^). Alternatively, one can model the GEV parameters using temporal or spatial 
covariates. A few examples include /i — fj,i + 112 • t, where t is time, or a — ai + a2 ■ 
lat -\- (73 • Ion + (T4 ■ elev, which considers effects of latitude, longitude, and elevation on the 
scale parameter. No matter the structure of the parameter 0, define the density function 
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gim; 0) = -£^G{m; 0). Then, the maximum hkehhood estimate of is 

^MLE = argmax^ Jj5f(m | 0) 



(4) 



This maximization is often done numerically, and has been implemented in a number of 
software programs including R (R Development Core Team, 2010) using the function fgev 
in the package evd. The density function for the fitted model is obtained by plugging in the 



2.2. Pricing a Contract Through Simulations 

Once we have a fitted model, we can use this to estimate the necessary pure premium 
and risk loading. This requires estimates of the first two moments of the unknown payment 
variable L, as 



where L{m]s,t)'^ is the loss payment for realization M = m raised to the rf*^ power {d = 1 
or 2), and g{m) is the density function of the maxima. The first type of weather derivative 
discussed has L(m; s, t^ — a'^ for m> s, which means the integral can be evaluated exactly 
as 



The second and third types of derivatives involve more complicated integrals, so we use 
monte carlo techniques to estimate them. This approach can be used to estimate moments 
for other types of derivatives with even more complicated payment structures, and is thus 
the most general approach. Here we draw a large iid sample Mj ~ G{m) for i — 1, /, and 
for each draw we compute the payment L{mi). Assuming that £'(|L(M)p) is finite, then by 
the Strong Law of Large Numbers as / — > oo, sample means converge to the first and second 
moments as (Robert, 2007) 



maximum likelihood estimate as g{m;(f)MLE)- 





(5) 




i=l 



(almost surely) 



(6) 
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and 

1 ^ f 

-^L{mif ^ E{L{Mf) = I L{mfg{m)dm (almost surely). (7) 

Furthermore, if the fourth moment is finite, as J L{mYg{m) dm < oo, then by the Central 
Limit Theorem we know that the sample average in equation 6 is asymptotically normal 
with variance var(L(M))//, and the sample average in equation 7 is asymptotically normal 
with variance var(L(M)^)//. We estimate the expected payments under the second and 
third contracts by drawing Mi,...,M/ ~ G{m \ 4>), and compute L{rni) and L{rniY f*^^ 
each of draw using / = 1,000,000 total draws. This total number was chosen to provide 
highly accurate estimates within a reasonable time, and simulations showed it was sufficiently 
high to eliminate concern for purely numerical monte carlo error. Sample averages of each 
converge to the theoretical first and second moments, which are used in the pricing model. 
One example of a risk-loaded premium based on marginal variance is 

1 ^ 

P = E{L) + R{L) = y E ^^'"^) + ^ ■ 

for some dollar amount A, chosen to satisfy some risk tolerance criteria. 

2.3. Example: Extreme temperature m Phoenix, Arizona 

As an example of how this model may be used, consider pricing a weather derivative 
with payments whenever the maximum daily summer temperature in the city of Phoenix, 
AZ exceeds some high threshold s. On June 26, 1990, Phoenix airport was forced to close be- 
cause the temperature exceeded 122 degrees Fahrenheit. Aircraft operating manuals did not 
provide information for takeoff and landing procedures in temperatures above 120 degrees 
Fahrenheit. The closure caused the predictable sort of economic disruption which accompa- 
nies airport closures. Wc envision a weather derivative as a useful tool in this situation. 

To price the derivative, we collect maximum daily summer temperatures at the Phoenix 
airport, for year i and day j, where j — 1, 92 (the 92 days in June, July, and August) 



i=l \ i=l 



mi) 



(8) 
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Maximum Daily Summer Temperature Autocorrelation 




1940 1960 1980 2000 5 10 15 

Vfear Lag 



Figure 1: Left: Maximum annual summer temperature at Phoenix International Airport from 1933-2010 
(with some missing values). The dashed line shows the annual trend, with statistically significant positive 
slope (p-val = 0.007). Right: Empirical autocorrelation of maximum daily summer temperatures. There is 
no evidence annual maximum daily temperatures are autocorrelated, as the value at all lags greater than 1 
falls below the 95% confidence interval line obtained from white noise sequences. 

for years 1933 to 2010. This data comes from the National Chmate Data Center. For each 
year i, we take the block maximum rrii = max(?/j i, 2/1,92), and model these annual maxima 
Generalized Extreme Value distribution. Plotting these data, we observe evidence of 
a slight positive trend over time (figure 1). A simple linear model of maximum temperature 
versus year shows a statistically significant positive slope of 0.03363, with p-value 0.007. We 
also find no evidence that annual maximum temperatures are autocorrelated. 

We estimate the GEV parameter (p using maximum likelihood estimation, as shown in 
equation 4, but with the possibility of a trend on the location parameter, as /i = /xi + /i2 ■ t 
where t is year. Thus, the GEV parameter here is actually = (/ii, /i2, o", 0- The maximum 
likelihood estimates (with standard errors shown in brackets) are /ti = 113.367 [0.250], 
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Table 1: Phoenix, AZ example: Model first and second moments for the type 1 weather derivative with flat 
payment L = 1000 paid whenever the maximum daily temperature M > s in the year 2011, using equation 5. 



Threshold s 


114 


116 


118 


120 


122 


124 


E{L) 
E{L'') ■ 10^3 


759.11 
759.11 


391.84 
391.84 


144.11 

144.11 


41.61 
41.61 


9.72 
9.72 


1.79 
1.79 



Table 2: Phoenix, AZ example: Model first and second moments for the type 2 weather derivative with 
proportional payment L = 1000 • (m — s) in the year 2011, for varying thresholds of s. Estimates are based 
on 7 = 1, 000, 000 monte carlo draws used with using equations 6 and 7. 



Threshold s 


114 


116 


118 


120 


122 


124 


E{L) 
E{L^) ■ 10-^ 


1,882.13 
7,336.56 


732.20 
2,369.34 


224.57 
627.82 


56.39 
137.98 


11.59 
24.45 


1.87 
3.26 



/(2 = 0.035 [0.011], 0- = 1.931 [0.176], and I = -0.090 [0.078]. Figure 2 shows some common 
diagnostics and the return level plot. The return level plot shows the expected number of 
years before an exceedance of a certain level is reached. This is the same as the reciprocal 
of the probabiUty of a specified exceedance, and forms the basis for statements such as 
describing an event as a "once every 50 years" event. 

With the fitted model for maximum summer temperature, we can estimate the first 
and second moments of various weather derivative payments in the year 2011. Estimated 
moments for the three types of derivatives from section 1 are shown in tables 1, 2, 3, and 4. 
As the limit i — > oo, payments under the second and third types are equal. 

Tables hke these can be used to price a wide range of weather derivatives. Consider a 
weather derivative with payment 1000- (M- 118) for M < 125 and 7000 for M > 125, where 
M is the maximum summer temperature in Phoenix. Using equation 8 with A = 0.0001, the 
tables show the pure premium should be 223.89 + 0.0001 • (618.16 ■ 10^ - 223.89^) = 280.69. 
Premiums for other limits, strike values, and payment structures can be estimated from the 
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Probability Plot 



Quantile Plot 
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Figure 2: Diagnostics from the maximum likelihood fit to the Phoenix summer temperature data. Top left: 
comparison of empirical and model probabilities. Top right: comparison of empirical and model quantiles. 
Bottom left: The return period is the expected number of years required for the process to exceed the 
corresponding return level. Bottom right: Model density function for 2011 maximum summer temperature 
in Phoenix AZ. 



12 



Table 3: Phoenix, AZ example: Model first moments for the type 3 weather derivative with proportional 
payment L = 1000 • (m — s) up to limit 1000 ■ {t — s) in the year 2011, for varying thresholds of s and t. 
Estimates are based on / 1, 000, 000 monte carlo draws used with using equations 6 and 7. 



t \ 


114 


116 


118 


120 


122 


124 


119 


1,766.86 


616.93 


109.30 








121 


1,855.93 


705.99 


198.37 


30.18 






123 


1,877.34 


727.41 


219.78 


51.59 


6.80 




125 


1,881.44 


731.51 


223.89 


55.70 


10.90 


1.18 


oo 


1,882.13 


732.20 


224.57 


56.39 


11.59 


1.87 



Table 4: Phoenix, AZ example: Model second moments for the type 3 weather derivative with proportional 
payment L = 1000 • (m — ,s) up to limit 1000 ■ {t — s) in the year 2011, for varying thresholds of s and t. 
Estimates are based on / = 1, 000, 000 monte carlo draws used with using equations 6 and 7. Values shown 
have order 10^. 



t 


114 


116 


118 


120 


122 


124 


119 


5,894.14 


1,383.33 


98.21 








121 


6,914.34 


2,050.63 


412.61 


26.27 






123 


7,243.24 


2,294.71 


571.86 


100.69 


5.84 




125 


7,321.98 


2,357.23 


618.16 


130.77 


19.70 


0.97 


oo 


7,336.56 


2,369.34 


627.82 


137.98 


24.45 


3.26 
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same general approach once a fitted model has been obtained. 

3. Background on Spatial Statistics and Spatial Extremes 

To extend beyond a single location, we need to consider models for multivariate and 
spatial extremes. Copulas provide a useful tool for modeling joint dependence, but they often 
fail to model extremes well (Mikosch, 2006). Since weather has a natural spatial domain, a 
better choice is to build spatial models designed for extremes, and use those to determine 
the joint dependence in payments in a collection of weather derivatives. In this section we 
present some background on spatial statistics, max-stable processes, and statistical methods 
for fitting max-stablc processes to data. The goal is to convey the benefits and motivation 
for using max-stable processes, and to outline the statistical method for fitting max-stable 
processes to data implemented in the R package SpatialExtremes. Readers interested in 
more details of spatial statistics and max-stable processes can find much greater explanation 
in Cressie (1993), Schlather (2002), and Padoan et. al. (2010). 

3.1. Background on Spatial Statistics 

The basic object in spatial statistics is a stochastic process Y{x),x e X where X is a 
subset of BP, usually with p — 2. Let 

S{x) = E{Y{x)), xeX 

be the mean of the process defined for all of X, and assume that the variance of Y{x) 
exists everywhere in X. The process is said to be Gaussian if for any K > 1 and locations 
xi, ...,xk, the vector {Y{xi), Y{xk)) has a multivariate normal distribution. The process 
is strictly stationary if the joint distribution of {Y{xi), ...,Y{xk)) is the same as {Y{xi + 
h), ...,Y{xk + h)) for any h & X and for any K points xi, ...,xk- For a Gaussian process, 
strict stationarity implies 

Coy{Y{xi),Y{x2)) — C{xi — X2) for all xi,X2&X 
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That is, the covariance of the process at any two locations is some function C which depends 
only on the separation vector between points, and not the particular locations. This is also 
called second-order stationarity. Next, we define the variogram through the relation 



where the quantity 27 is the variogram, and 7 is the semi-variogram. Under the assumption 
of strict (or second-order) stationarity. 



where p{h) is the correlation between two locations separated by vector h. Further, if we 
have 7(/i) = 7(| |) for all /i e X, meaning if the semi-variogram only depends on h through 
its length \\h\\, then the process is isotropic. The correlation function p{h) is then usually 
chosen from one of the valid famihes of correlations for Gaussian processes. A few common 

choices are Whittle-Matern, 



where Ci, C2 and u are the nugget, range, and smooth parameters, F is the gamma function 
and Ki, is the modified Bessel function of the third kind with order u. It is common to fix 
the nugget as ci = 1, which forces p{h) — )■ ci = 1 as ^ 0. This is a reasonable assumption 
for many environmental processes, and we make this assumption throughout this paper and 
do not attempt to model the nugget. Throughout the remainder of this paper, the unknown 
spatial dependence parameter is called 9 — (c2, i^). 



Var(y(xi) - Y{x2)) = 27(^1 - X2) 



^{h)^C{0)-Cih) = Cmi-pih)) 




Cauchy, 




and powered exponential 
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3.2. Multivariate Extreme Value Distribution 

The final definition we need before we can introduce the model for spatial extremes is 
the multivariate extreme value distribution. Let (Fji, Fj^-), i = l,...,n be independent 
and identically distributed replicates of a fT— dimensional random vector and let M„ = 
(M„i, Mnx) be the vector of componentwise maxima, where M„jt — max(Yifc, Ynk) for 
k — 1, K. A non-degenerate limit for exists if there exist sequences a^^ > and bnk, 
k — 1, K such that 

f M^i - Ki ^ -hnK ^ \ ^, ^ 
hm F I < mi, < mx = (^{'mi, rriK)- 

Then G is a multivariate extreme value distribution, and is max-stable in the sense that for 
any n > 1 there exist sequences A^k > 0, Bnk, k = 1, ...,K such that 

The marginal distributions of a multivariate extreme value distribution are necessarily uni- 
variate Generalized Extreme Value distributions. 

3.3. Max-stable Processes 

Here we introduce the spatial analog of the multivariate extreme value distribution. Let 
Z{x),x G X C RP he a stochastic process. If for all n > 1, there exist sequences a„(a;), bn{x) 
for some xi, xk G X such that 

p f Z{xk)-bn{xk) ^ ^^^^^^^ ^ ^ ^\ ^ {,^x,), z{xk)) 

n-^oo \ an{Xk) ) 

then G'a;^^...^a;^ is a multivariate extreme value distribution. If the above holds for all possible 
xi, Xk £ X for any if > 1, then the process is a max-stable process. To briefiy summarize, 
if we have a max-stable process Z{x) defined for all x e X, then at any single location x ^ X 
the distribution of Zix) is GEV, and all finite vectors (Z(x})^ Z{xk)) follow a multivariate 
extreme value distribution. In this sense, max-stable processes arc the infinite dimensional 
generalization of multivariate extremes, and nicely extend the GEV to spatial domains. 
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Constructing a max-stable process is accomplished through a point process approach. 
Let Y{x) be a non-negative stationary process on Rp such that E{Y{x)) = 1 at each x. Let 
n be a Poisson process on i?+ with intensity dw/w"^. If Yi{x) are independent rephcates of 
Y{x), then 

Z[x) = maxWiYi^x), x G X 

i 

is a stationary max-stable process with unit Prechet margins (De Haan, 1984). Prom this, 
the joint distribution may be represented as 

P{Z{x) < z{x),x e X) =ex-p\-E ( sup ) j (9) 

L \xex z{x) J ) 

Varying the choice of the process Y(x) gives different max-stable processes. Smith (un- 
published manuscript, 1990) constructed a process known as the Gaussian extreme value 
process by taking Yi{x) to be a multivariate Gaussian centered at the point Xi with covari- 
ance matrix S. Smith also introduced the "rainfall-storms" interpretation in 2 dimensions: 
think of as the space of storm centers, Sj as the magnitude of the i*^ storm, and Yi{x) 
as the shape of the storm centered at position Xj. The maximum of independent storms at 
each location x is taken to be the max-stable process. A realization of this process is shown 
in figure 3. A particular strength of the Smith model is the ability to handle anisotropy, as 
shown in the figure. This comes from the off diagonal covariance parameter E12 in 




S12 E22 



We will consider the use of the Smith model only to check the assumption of isotropy, which 
is required for the next class of max-stable processes. Schlather (2002) introduced a flexible 
set of models for max-stable processes, termed extremal Gaussian processes. Consider a 
stationary Gaussian process Y{x) on R^ with correlation function p{-;9) and finite mean 
S — £'max(0, Y{x)) E (0, 00). Let Wi be a Poisson process on (0, 00) with intensity measure 
S~^w~'^ds. Then 

Z{x) = max Wj max(0, 1^(3;)) 

i 
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is a stationary max-stable process with unit-Prechet margins. The bivariate distribution 
function is 

1/2- 



P(Zi < zi,Z2 < Z2) = exp 





"1 ^ r 




H 


{-I 















l + a-2(p(/i;^) + l) 



Z1Z2 



{zi + Z2y 



(10) 



where p{h; 9) is the correlation of the underlying Gaussian process Y{x) and h — \ \xi — X2\\- 
Figure 3 shows one reahzation of a process with the Whittle-Matern correlation function. 
These processes are flexible, and produce plausible realizations of environmental processes. 
We have elected to concentrate on the Schlather model in this paper. 



Smith Model 



Schlather Model 




2 4 6 8 10 




Figure 3: Left: A realization of a Smith process with strong anisotropy. Right: A reahzation of a Schlather 
process with Whittle-Matern correlation and parameter = [ci = 1, C2 = 3, = 1). 



3.4- Maximum Composite Likelihood Estimation for Max-stable Processes 

A potential stumbling block to using max-stable processes is that the closed-form expres- 
sion of the joint likelihood shown in equation 9 can only be written out for dimension one 
or two. This means only the univariate likelihood (which is GEV) and bivariate likelihood 
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function (as shown in equation 10) are available in closed form. When one considers the joint 
distribution of a max-stable process at three or more locations, a closed-form expression for 
the joint likelihood is unavailable. One way to proceed with a likelihood-based approach is 
to substitute the composite likelihood for the unavailable joint likelihood. The composite 
log-likelihood is defined as 

TV 7-1 I 

ecie;z) = EE E ^Ogfiz,,n,Zj,n;e) (11) 
n=l i=l j=i+l 

where each term f{zi^n, Zj^n', is a bivariate marginal density function based on locations i 
and j. The two inner sums sum over all unique pairs, while the outer sums over the N i.i.d. 
replicates. Similar to the full likelihood function, the parameter which maximizes a com- 
posite log likelihood can be found, and is termed a maximum composite likelihood estimate, 
or MCLE. The maximum composite hkelihood estimator is consistent and asymptotically 
normal (Lindsay, 1988) (Cox and Reid, 2004) as 

OMCLE^N{eJ) with i ^H{e)j-\e)H{e). (12) 

where H{e) = E{-Helc{0] Z)) is the expected information matrix, J{e) = V{De£c{0] Z)) is 
the covariance of the score, and Hq is the Hessian matrix, Dq is the gradient vector, and V is 
the covariance matrix. In the setting with the full likehhood and MLE, we have H{9) — J (6), 
but for the composite likelihood setting these matrices are not equal. 

Padoan et. al. (2010) used composite hkehhoods to model the joint spatial dependence of 
extremes, and implemented their work in the R package SpatialExtremes. The maximum 
composite likelihood estimator Omcle is found using the numerical optimizer command 
optim. Estimates of the variances are found by plugging in 9mcle into expressions for H 
and J: 

N I-l I 

H{9mCLe) = - E E E log /(^n,i, Znj] OmCLe) 
n=l i=l j=i+l 
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N I-l I 

A^MCLe) = "X^X^ X] Dolog f{Zn,i, ZnJ, OMCLE)D0log f{Zn,i, ZnJ, OmCLe)'^ 
n=l 1=1 j=i+l 

Model selection is based on minimizing the composite likelihood information criteria 
(CLIC) (Varin and Vidoni, 2005), equal to 

—2ic{9MCLE', Z) — ti (^J{9mcle)H{9mcle) (13) 

where the second term is the composite log-likelihood penalty term. 

In our setting, this approach is used as follows: we begin with Y independent and 
identically distributed reahzations of an observed set of spatial extremes data, for loca- 
tions Xi,...,Xk- Thus there arc Y replicates and K locations. For each location x^, we 
transform the GEV data to unit-Frechet margins by first estimating all GEV parameters 
fi{xk)i o'{xk),i{xk), k = 1, ...,K, then using these in the transformation shown in equation 3. 
Next we obtain the maximum composite likelihood estimate for the dependence parameter 
9mcle — (c2, i^) of the max-stable process using the composite likelihood approach outlined 
above. The result is a fitted model for the extremes at the K specific locations, with spatial 
GEV parameter = (/i(xjfc), a"(xfe), ^(x^). A; = 1,...,K) and spatial dependence parameter 

4. Simulating Losses from Extremes at Multiple Locations 

Prom the fitted model, we can simulate a collection of max-stable process Zi{xk),i — 
for the same k — 1,...,K locations as the observed data, and then transform back 
to the original scale of extremes data by transforming margins using the estimated GEV 
parameter 0. From this we can compute the payments from weather derivatives at each 
location as L{mi^k] s, t). 

To price a collection of K weather derivatives, with jointly dependent losses, there are 
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several quantities of interest. First, the total variability of loss payments is 

(K \ K K-l K 

^Lk \ + 5Z 2 •cov(Lfe,Lfc/). (14) 

k=l J k=l k=l k'=k+l 

Now consider a portfolio of X — 1 derivatives, with the seller deciding whether or not to 
write a K*'^ derivative. The additional derivative will increase the total portfolio variance by 

MVk = var (j^ Lkj - var (j^ Lk] = var(L^) + J] 2 • cov(Lfe, Lk). (15) 
\fe=i / \fe=i / k=i 

The covariance for any two derivatives at locations Xk and Xk' is 



cov(Lfc, Lk') = E{Lk ■ Lk') - E{Lk)E{Lk') 



(16) 



Each of these quantities may be estimated from a large collection of / simulations. The 
total variance of the portfolio is 

f -I i A' / -I i A' \ ^ 



var 



.fe=i 



=1 k=l 



1=1 k=l 



(17) 



the marginal variance for adding a K^'^ derivative is 



K 



''K-l 

var I ^ L{mi^k) 
.fe=i 



K-l 



m3i{Lk) + XI 2 ■ cov(Lfc, Lk), (18) 



k=l 



and the covariance between any two derivatives is 



cov(Lfe,Lfc/)= ^Iv(mi,fc) • L(mi,fc/)j - ^y^L(mi,fc)j ^y^L(mi,fc/)j (19) 

4..1. Simulated Example 

We evaluate the approach through simulations, first with a single detailed case and then 
larger numbers of simulations. It will be convenient to define some new notation to keep simu- 
lation results clear. Call the true full parameter — {ii{xi), a{xi),^{xi), ij,{xk), cr{xK),^{xK), C2, i^), 
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and the estimated parameter 9. Call the true marginal variance MV{9), and call an estimate 
of this based on a fitted model MV{9). We begin with a single detailed case. 

We simulated a max-stable process with parameters chosen to mimic annual temperature 
maxima in North America. The process had unit-Frechet margins and Whittle-Matern 
covariance with dependence parameter icx^c^^v) = (1,3,1) for 75 years at 20 locations 
randomly placed on a 10 by 10 grid. Call the vertical dimension latitude (lat) and the 
horizontal longitude {Ion). To make this data consistent with annual temperature maxima, at 
each location we transformed to the GEV scale by specifying parameters //(x) = 110 — lat /2, 
a{x) — 1.5 + lat/5, and ^{x) — —0.1. The basic idea was to imagine higher latitude locations 
having overall lower extreme temperatures, but higher variability of extremes. We used 
these transformations for each of the 20 locations to produce a max-stable process with 
GEV{/j{x),a{x),^{x)) margins. We fix this as the "observed" data. 

Next, we analyzed these data using composite likelihood estimation. For each location 
Xk, we obtained a maximum likelihood estimate 4>{^k) = f'^{xk),o-{xk),i{xk),k = 1,...,K 
using equation 4, and then used these to transform each margin to unit-Frechet using equa- 
tion 3. We fit a max-stable process with Whittle-Matern correlation with nugget parameter 
1, and obtained maximum composite likelihood estimate (£2, i>). Using this fitted model, we 
simulated a large number of processes and transformed them to the temperature scale using 
jl{xk),o-{xk),i{xk) at each location Xk- 

Thus we have described a means of simulating i = 1,...,/ extreme temperature events 
rrii^k for locations Xi,...,Xk from our fitted model. This information was used to estimate 
payments for weather derivatives by computing L{miy, s,t) . Table 5 shows results from 
one simulation. We simulated 200,000 extreme temperature events at the same 20 locations, 
and used these to compute payments Li^k for ^ = 1, 200, 000 and k — 1,...,K. Shown 
are payments for a weather derivative paying 1 when T > 112, the aggregate payment 
^^k=i s-nd the payments for a possible weather derivative L2Q. The marginal variance of 
adding the 20*'* derivative was estimated as MV2q{9) = 22.788 - 21.216 = 1.572, which is 
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Table 5: Simulated payments for weather derivatives paying 1 when T > 112 using the fitted model with 
parameter 0. This information shows the marginal variance for adding the 20*^* policy is MVw{0) = 1.572. 
Using the true model, MV2o(^) is 1.250. The difference comes from parameter estimation error in 9. 



Event 


Li 


L2 . 


Ll9 


J2f=i Lk 


L20 


J2kLi 


1 




















2 


1 





1 


8 


1 


9 


3 


1 


1 





4 





4 


200,000 











1 





1 


Mean 


0.095 


0.130 . 


. 0.189 


2.907 


0.089 


2.996 


Variance 


0.086 


0.113 . 


. 0.153 


21.216 


0.081 


22.788 



clearly much larger than 0.081, the estimated variance of the payments when dependence 
terms are ignored. An additional 200,000 simulations from the true model with parameter 

9 shows the true marginal variance MV2o{9) is 1.250. This particular simulation showed an 
overestimation of the marginal variance of 25.7%, which is clearly a substantial error, but 
not when compared to the error from ignoring spatial dependence. 

4.2. Simulation Study of Performance 

We evaluated the performance of this method in estimating the marginal variance of 
adding a 4*^* weather derivative to an existing portfolio composed of Li, L2. and L3. This 
quantity is key to pricing a risk load for L4. We randomly placed K — 25 locations on a 10 by 

10 grid, and randomly selected 4 of these to represent locations of weather derivatives. The 
target quantity was MV4, the marginal variance of adding a fourth derivative. We estimated 
this quantity using two methods: 

1. Estimate MV4 using equation 18, which accounts for spatial dependence by fitting a 
max-stable process and uses simulations from the model, with fitted parameter 6 — 
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(//(xi), a{xi),i{xi), ji{x4), a{x4),i{x4), 62, i>). 
2. Estimate MV^ using equation 17, which fits a GEV to the data at location k — A but 
does not account for spatial dependence among the derivatives, with fitted parameter 

is ^= (/i4, 0-4,0 • 

Again, call the true full parameter 6 — a"(xi), ^(xi), //(x^), <7(a;x), ^(^^-ftr)) C2, i^), 

and the estimated parameter 9. Call the true marginal variance MV{9), and an estimate 
MV{9). The true marginal variance was found by simulating I — 1, 000, 000 reahzations of 
a max-stable process under the true parameter 9, and using equation 18. Method 1 uses the 
same approach, but with estimated parameter 9, as we showed in the single example above. 
Method 2 ignores spatial dependence. The first measure of error we use is percentage error, 

MV,{0) - MV0} 



where J = 1, ...,500 refers to a simulation run. This choice preserves the sign of estimation 
error. Results are shown in figure 4. Here, we see the peril of ignoring spatial dependence 
of losses in a collection of weather derivatives. The right column shows that as the range of 
the spatial dependence increases, the underestimation bias of estimating marginal variance 
MVi increases. The left column shows the unbiased results obtained from incorporating 
dependence using the method of this paper. 

We also show the asymptotic results for Method 1 in table 6. Here, we use a slight variant 
of estimation error called mean absolute percentage error, 

I J^\MVA9) - MVA9)\ 
MAPE^-Y : ''Jy^^g/'" - (21) 

This choice does not preserve the sign of error, but is more suited to showing asymptotic 
results. Results from 150 simulations in a variety of years and dependence ranges are shown 
in table 6. For all dependence ranges shown, the error in estimation falls as more data is 
available. The remaining error is primarily due to parameter estimation error in 9 . 

24 



Method 1 : Range = 0.5 



Method 2: Range = 0.5 



-1.0 



-0.5 0.0 0.5 

Percentage Error 

Method 1 : Range = 3 



-0.5 0.0 0.5 

Percentage Error 

Method 1 : Range = 8 



-0.5 0.0 0.5 

Percentage Error 



-0.5 0.0 0.5 

Percentage Error 

Method 2: Range = 3 



-0.5 0.0 0.5 

Percentage Error 

Method 2: Range = 8 



-0.5 0.0 0.5 

Percentage Error 



Figure 4: Comparison of methods 1 and 2 for estimating the marginal variance MV4 under three spatial 
dependence scenarios, with 500 simulations in each. The heavy line at signifies the true MV4{9). The top 
row corresponds to a short-range dependence process, the middle row is medium-range, and the third row 
shows long-range spatial dependence. The left column uses the approach outlined in this manuscript, also 
called Method 1, which incorporates spatial dependence. The right column shows the non-spatial model in 
Method 2. Results are plotted as percentage deviation from the true marginal variance, using equation 20. 
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Table 6: Mean absolute percentage error {MAPE) in estimating MV4, the marginal variance of adding a 
fourth weather derivative to a portfolio (using only Method 1, which incorporates the spatial dependence). 
Estimates are averages from 150 simulations each based on 25 locations, and are computed using equation 21 
for 50, 100, 250, and 500 years of data. For each spatial dependence range shown, error falls as the number 
of years of data increases. The remaining error can be attributed to parameter risk. 



Range 


Y=50 


Y=100 


Y=250 


Y=500 


Short 0.5 


0.160 


0.121 


0.071 


0.051 


Medium 3 


0.233 


0.139 


0.108 


0.064 


Long 8 


0.231 


0.151 


0.087 


0.064 



5. Using Simulation Output to Price Weather Derivatives 

5.1. Incorporating Marginal Variance into Risk Premiums 

Here we discuss using the models to price a risk load. Feldblum (1990) described five 
methods for an insurer to determine the risk load for writing a policy with unknown loss 
L. The first two methods discussed were the variance approach, where R = X ■ var(L), and 
the standard deviation approach, where R = \ ■ \J Kreps (1990) and also Philbrick 
(1991) discussed how a new policy adds risk through the marginal increase in total variance, 
shown in equation 18. Gogol (1992) cautioned that these methods of pricing risk loads are 
order- dependent if losses are correlated, leading to a mismatch between individual renewal 
risk loads and the total portfolio risk load, best illustrated by a toy example. 

Consider two correlated policies, Li and L2, and consider computing the risk load based 
on marginal variance in two ways: 

1. A risk load using marginal variance for the total loss L — Li + L2 would be 
A (var(Li) + var(L2) + 2 • cov(Li, L2)). 

2. A risk load computed individually would go as follows: when L2 is priced, the risk load 
is A(var(L2) + 2 • cov(Li,L2)). When Li renews, it receives risk load A(var(Li) + 2 • 
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cov(Li, I/2)). The sum of these renewal risk loads is A (var(Li) + var(L2) + 4 • cov(Li, L2)), 
which has double counted the covariance terms and does not match the total portfolio 
risk load. 

The example demonstrates the danger in careless accounting of covariance terms. The ap- 
proach we take to pricing a portfolio of dependent weather derivatives follows the work of 
Mango (1998). In Mango's terminology, we use the covariance-share method, which appor- 
tions the total covariance between policies Lj and and computes risk loads as 

R{Lk) = A ^var(Li^) + 2 ^ aj,K ■ cov(L^, L,) j (22) 

for any < a^ x < 1. This quantities a^ x are chosen to split the respective covariance 
terms and ensure the sum of individual renewal risk loads matches the total portfolio risk 
load. One reasonable choice splits the total covariance in proportion to the expected losses 
of policies j and K, as 

- E{L,) + E{L,y 

Under this choice, we always have aj^x + clkj = 1, so the risk loads will always be renewal- 
additive. Relevant quantities are estimated from large numbers of event simulations, and 
the risk load is 

R{Lk) = A (vSxiLK) + 2 ^ aj^K ■ cov(Lk, L,) j (24) 

using equations 17 and 19, where 

. _ 1 E?=i L{rni,K) , . 

I Ei=i + I Ei=i Hm,K) 

5.2. Example: Midwest Temperature Data 

We illustrate the methodology on US temperature data. The data, freely available from 
the National Climate Data Center (http://cdiac.ornl.gov/ftp/ushcn_daily/) .come 
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from 39 locations in the midwestern United States with complete summer (June 1 - Au- 
gust 31) temperature records from 1895 to 2009. All sites are located between -93 and -103 
degrees longitude, and 37 to 45 degrees latitude, shown in figure 5. We use all 39 locations 

to estimate the max-stablc process, but we only consider weather derivatives at 4 of these 
locations, labeled 1-4 and drawn with triangles on the figure. Call the maximum summer 
temperature at these k = 4 locations Mj^^, with payments Li^k defined as 

1. Li,i = 1000 if {Mi,i > 107} and otherwise 

2. Li^2 = 300 • (Mi,2 - 105) when {105 < M2 < 110}, 1500 when {Mi,2 > 110}, and 
otherwise 

3. Li,3 = 200 if {Mi,3 > 105} and otherwise 

4. Li,4 = 200 if {Mi,4 > 102} and otherwise 

The application proceeded in two steps. The first is to use data from all 39 locations to 
fit a max-stable process in the study region, and the second is to then simulate temperature 
events from the fitted model only at locations 1-4 to estimate the renewal-additive risk load 
and premium for adding a weather derivative at location 4. 

To investigate the possibility of a trend in maximum daily temperatures over time, we 
fit simple linear models to maximum daily temperature versus year, but found only 4 out of 
39 locations showed statistically significant slopes at the p — 0.01 level (this lower level was 
selected to reduce the false-positive rate which occurs with multiple tests). Furthermore, all 
four slopes were negative. We also fit GEV models to data from each station allowing for 
a time-varying location parameter as Hk = A*fc,o + A*fc,i • t, where t is year, but found only 7 
differed significantly from (again at the p=0.01 level), and again, all were negative. These 
locations were spread throughout the study region, and showed no discernible spatial pattern 
or clustering. We concluded that there was no evidence of a widespread shift in maximum 
temperatures over time throughout the entire region, and dropped the time-varying GEV 
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location parameter. However, just as a precaution we also conducted a separate analysis of 
the data including these 7 negative trends, but found it had little impact on the results. 

We fit ordinary GEV models to each station, and obtained maximum likelihood estimate 
= {fi{xk),o'{xk),i^{xk)) for k = 1,...,39. Diagnostics like those shown in figure 2 gave no 
indication the GEV was inappropriate for any of these locations. These fitted models were 
used to transform data at each location to unit-Frechet. Next we assessed the appropriateness 
of using a max-stable process for the dependence. We first fit a Smith process to the unit- 
Prechet data to check for anisotropy, but did not see strong evidence of anisotropy. The 
parameter estimate of covariance E were En = 2.064 [0.020] E22 = 1.897 [0.020], and 
^12 = —0.085 [0.009] 0, where the number in brackets is the standard error of the estimate 
(when Ell = E22 and E12 = 0, we have perfect isotropy). We next considered the more 
flexible Schlather model with Whittle-Matern, Cauchy, and powered exponential correlation 
functions, and found the Whittle-Matern to be the best with the lowest CLIC score. Using 
the Whittle-Matern correlation model, we obtained maximum composite likelihood estimates 
of the range and smooth as C2 = 4.6819 [1.2975] and z> = 0.3155 [0.04625], where the number 
in brackets is the standard error of the estimate. 

Next, we simulated / = 100, 000 max-stable processes from our fitted model at the 
four locations with weather derivatives. Using the GEV estimates {fi{xk),a{xk),^{xk)) for 
/c = 1,2,3,4 we transformed the unit-Prechet margins to GEV at each location. Thus, we 
had simulations of maximum summer temperatures Mi^k for i = 1, 100,000 at the k = 4 
locations. From these, we computed the payments Lj^^ under the four contracts considered. 
Table 7 shows a few of these simulations. 

From equation 24 and the information in the table, we compute the risk load for contract 
L4 as 

R{Li) = (46.95 + 2 (0.8229 • 8.43 + 0.3636 • 29.93 + 0.1995 • 28.46)) • 10^ • A = 93, 044 • A. 
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Table 7: Payments for weather derivatives in the Midwestern temperature example. / — 100, 000 simulated 
extreme temperature events are simulated at locations 1-4 in figure 5. The covariance share quantities aj^K 
are estimated using equation 25. 



Event 


Li 


L2 


Ls 




L4 




1 




















2 





757.76 





757.76 





757.76 


3 




















4 


1000 


964.02 





1964.02 


444.94 


2408.96 


100,000 




















Mean 


221.75 


96.751 


11.892 


330.393 


55.271 


385.664 


Variance (-10"^) 


172.58 


99.89 


6.35 


381.38 


46.95 


561.98 


Cov(Lfe,L4)(-10-3) 


28.46 


29.93 


8.43 


66.82 








0.1995 


0.3636 


0.8229 
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Midwest Temperature Example Locations 




-102 -100 -98 -96 
Longitude (degrees) 



Figure 5: Locations of the 39 stations in the Midwest temperature example used to fit the max-stable 
process to maximum summer temperature. The locations in triangles labeled 1-4 are the places where 
weather derivatives are priced. 



This is roughly half of the total increase of (561.98 - 381.38) • 10^ • A = 180, 600 • A. The 
remainder would be apportioned to the risk loads of first three derivatives as they renew. 

When we included the 7 time-varying location parameters, we computed a risk load of 
90, 914 • A, a reduction of only 2.3%. This alleviates concerns that we might have wrongly 
ignoring trends in the GEV location parameter fj,. If the inclusion of trends on the GEV 
location parameters fik,k = 1,...,39 had resulted in a substantially larger risk estimate, it 
might warrant the inclusion of trends as the more conservative choice. However, the limited 
statistical evidence of trends combined with such a small reduction in the risk estimate 
supports dropping them altogether. 

6. Discussion 

We have described a means of pricing a collection of extreme weather derivatives based 
on simulations from max-stable processes. Naturally, there will be some error between the 
collection of simulated payments and actual payments. We discuss the errors introduced 
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from model selection, simulations, and parameter estimation. 

We have taken the approach to modeling spatial extremes as max-stable processes with 
Generalized Extreme Value margins, and naturally this model may not be appropriate for 
some spatial extremes data. For weather derivatives with payments based on maxima (or 
minima) of some weather variable, models based on block maxima (or minima) of the data 
make the most sense, and certainly the GEV distribution has appealing asymptotic properties 
for these data. Coles (2001) discusses diagnostics to check the validity of the GEV for 
the marginal data. Max-stable processes very naturally extend the GEV to the spatial 
domain, and are thus the logical choice for spatial block maxima data. While a goal is to 
extend the approach presented in this manuscript to include non-stationary fields, at present 
this approach can only handle stationary fields. Within the class of stationary max-stable 
processes, there are some choices of models. One can model the GEV parameters fi, a, 
and C, with spatial, temporal, or other covariates, and one can consider different correlation 
functions p{h) for the spatial dependence of the max-stable process. In this paper we did 
not show much detail on model selection, however the paper by Padoan et. al. (2010) shows 
the use of composite likelihood information criteria to handle model selection questions like 
these. 

The computational cost of simulations from a max-stable process is minimal, and thus 
one can simulate hundreds of thousands or millions of events with relative ease. Errors 
arising from numerical approximation in estimating the moments of payments assuming 
some fitted model using equations (6) and (7) are thus likely to be quite small, and can be 
made arbitrarily smaller with greater numbers of simulated events. 

The largest source of error in this approach is likely to come from parameter risk - that 
is, the error in estimating the GEV and max-stable process parameters and 9. Reducing 
parameter risk is best handled through fitting the process to more weather data: more years 
of data, more locations of data, or ideally both. A point worth stressing is that the data used 
to fit the max-stable process can (and probably should) contain far more locations than the 
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portfolio of weather derivatives. By adding additional points of data to fit the process, one 
reduces the parameter risk associated with estimating 6, the spatial dependence parameter. 

Our analysis is really a two-step procedure: the first transforms GEV margins to unit- 
Frechet by obtaining parameter estimate = {(ji,{xi),(j{xi).^{xi)^ ...^ (i{xk).,(y{xk).i{xk)), 
and then in a second step we fit a max-stable process to the transformed data to obtain 
dependence parameter estimates 9 = {62, u). We should point out that a single step pro- 
cedure is possible, and is implemented in the package SpatialExtremes, but this has two 
drawbacks. The first is that the numeric optimization of the likelihood needs to maximize a 
high dimension parameter. In our example on Midwestern temperature data, the dimension 
would be 39- 3-1-2 = 119 (and even larger if we kept time- varying GEV location parameters). 
The dimension raises concerns that the numeric optimizer may converge to a local maxima, 
not the global one. A second drawback is that single-step maximization of a max-stable 
process can be a painfully slow process, requiring orders of magnitude more time than a 
two-step procedure. With these drawbacks in mind, the two-step procedure was selected. 

One final comment is the potential mismatch between past and future weather extremes, 
particularly in the context of climate change. One can model the location /i and scale a 
parameters of the GEV with time covariates to allow for the possibility of non-stationary 
maxima in time. It is much less common to model the shape parameter ^ as anything other 
than a fixed number. We have illustrated the use of time covariates for modeling the location 
parameter as /i = fii + fi2't in the Phoenix airport temperature example. We caution readers 
not to extrapolate models such as these too far into the future. 
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